#!/usr/bin/perl -w
## quick way to extract accessbility distribution information on the two epitypes in Huy's epitype program's regions files
## author: Yaping Liu  lyping1986@gmail.com
##2013-11-6

##usage: perl extract_acc_in_epitype_nomeseq.pl cg_ref_coord.bw epi_summary.txt region_file_dir

use strict;
use Getopt::Long;
use File::Basename;

my $ref_cg=$ARGV[0] or die;
my $epi_summary=$ARGV[1] or die;
my $regions_file_dir=$ARGV[2] or die;
my $key=$ARGV[3];
my $output_epi_H = $epi_summary.".$key.epi_H.txt";
my $output_epi_L = $epi_summary.".$key.epi_L.txt";
my $methy_output_epi_H = $epi_summary.".$key.epi_H.methy.txt";
my $methy_output_epi_L = $epi_summary.".$key.epi_L.methy.txt";

#my %cg_hash=();
#open(CG,"<$ref_cg")  or die "no ref_cg coord files: $ref_cg\n";
#while(<CG>){
#	chomp;
	
#}
#close(CG);

open(HIGH,">$output_epi_H") or die "can't write epitype_H file: $output_epi_H\n";
open(LOW,">$output_epi_L") or die "can't write epitype_L file: $output_epi_L\n";
open(HIGHM,">$methy_output_epi_H") or die "can't write epitype_H METHY file: $methy_output_epi_L\n";
open(LOWM,">$methy_output_epi_L") or die "can't write epitype_L METHY file: $methy_output_epi_L\n";
open(EP,"<$epi_summary") or die "no epi_summary files: $epi_summary\n";
while(<EP>){
	chomp;
	next if($_!~/^chr/);
	my @splitin=split "\t";
	my $region_file=$regions_file_dir."/$splitin[0].$splitin[1]:$splitin[2]";
	
	my @gch_coords=();
	my @hcg_coords=();
	my $epi_1=""; ##default is Prob1 column.
	my $epi_2="";

	open(RG,"<$region_file") or die "no region files: $region_file\n";
	while(my $region_line = <RG>){
		chomp($region_line);
		
		my @regions_features = split "\t",$region_line;
		
		##find HCG coordinates.
		if($region_line =~ /^coor/){
			my $chr=$splitin[0];
			for(my $i=1;$i<=$#regions_features;$i++){
				my $start=$regions_features[$i]-1;
				my $end=$regions_features[$i];
				my $wl = `bigWigSummary $ref_cg $chr $start $end 1 -type=coverage 2>/dev/null | wc -l`;
				chomp($wl);
				$wl =~ s/(\d+)\s+\S+/$1/;
				if($wl == 1){
					push(@hcg_coords,$i);
				}elsif($wl == 0){
					push(@gch_coords,$i);
					#print "gch_coords:\t$i\t".scalar(@gch_coords)."\n";
				}else{
					die "weird cg coordinate..\n";
				}
			}
		}
		
		##count epi_1 epi_2's methylation level
		if($region_line =~ /^epitype_1/){
				my @ep1=();
				foreach my $i(@hcg_coords){
					push(@ep1,$regions_features[$i])
				}
				$epi_1=&mean(@ep1);
				if($epi_1 < 0){
					last;
				}
		}elsif($region_line =~ /^epitype_2/){
				my @ep2=();
				foreach my $i(@hcg_coords){
					push(@ep2,$regions_features[$i])
				}
				$epi_2=&mean(@ep2);
				if($epi_2 < 0){
					last;
				}
		}
		##count epi_1 epi_2's methylation level
		if($regions_features[0] eq $key){
			my @methys=split "",$regions_features[1];
			##output accessbility level
			my @accs=();
			foreach my $i(@gch_coords){
					#$i--;
					push(@accs,$methys[($i-1)])
				}
			my $acc_mean=&mean_binary(@accs);
			if($acc_mean >= 0){
				if($epi_1 > $epi_2){
					if($regions_features[2] > $regions_features[3]){
						print HIGH "$acc_mean\n";
					}else{
						print LOW "$acc_mean\n";
					}
				}else{
					if($regions_features[2] > $regions_features[3]){
						print LOW "$acc_mean\n";
					}else{
						print HIGH "$acc_mean\n";
					}
				}
			}
			##output methylation level
			my @ms=();
			foreach my $i(@hcg_coords){
					#$i--;
					push(@ms,$methys[($i-1)])
				}
			my $mc_mean=&mean_binary(@ms);
			if($mc_mean >= 0){
				if($epi_1 > $epi_2){
					if($regions_features[2] > $regions_features[3]){
						print HIGHM "$mc_mean\n";
					}else{
						print LOWM "$mc_mean\n";
					}
				}else{
					if($regions_features[2] > $regions_features[3]){
						print LOWM "$mc_mean\n";
					}else{
						print HIGHM "$mc_mean\n";
					}
				}
			}
			
		}

	}
	close(RG);
}
close(EP);
close(HIGH);
close(LOW);
close(HIGHM);
close(LOWM);

sub mean{
	my @ds=@_;
	my $sum=0;
	my $len=scalar(@ds);
	foreach my $d(@ds){
		$sum += $d;
	}
	if($sum <= 0){
		return -1;
	}else{
		return $sum/$len;
	}
	
}

sub mean_binary{
	my @ds=@_;
	my $methy=0;
	my $sum=0;
	if(scalar(@ds) <= 0){
		return -1;
	}
	for(my $i=0;$i<=$#ds;$i++){ 
		if($ds[$i] eq 'm'){
			$sum++;
			$methy++;
		}elsif($ds[$i] eq 'u'){
			$sum++;
		}
		
	}
	if($sum <= 0){
		return -1;
	}else{
		return $methy/$sum;
	}	
}
